function [N,Nxi,Neta,N2xi,N2eta,N2xieta,N2etaxi] = ShapeFunc(elem,p,z) 
% [N,Nxi,Neta] = ShapeFunc(elem,z)  
% N, Nxi, Neta: matrices storing the values of the shape functions on the Gauss points
%               of the reference element
%               Each row concerns to a Gauss point
% z: coordinates of Gauss points in the reference element
% elem: type of element (0: quadrilateral, 1: triangles)
% p: interpolation degree

xi = z(:,1); eta = z(:,2); 
if elem == 0
    if p == 1
        N    = [(1-xi).*(1-eta)/4, (1+xi).*(1-eta)/4, (1+xi).*(1+eta)/4, (1-xi).*(1+eta)/4]; 
        Nxi  = [(eta-1)/4, (1-eta)/4, (1+eta)/4, -(1+eta)/4]; 
        Neta = [(xi-1)/4, -(1+xi)/4,   (1+xi)/4,  (1-xi)/4 ]; 
        N2xi    = [zeros(size(xi)), zeros(size(xi)), zeros(size(xi)), zeros(size(xi))];
        N2eta   = [zeros(size(eta)), zeros(size(eta)), zeros(size(eta)), zeros(size(eta))];
        N2xieta = [1/4*ones(size(xi)), -1/4*ones(size(xi)), 1/4*ones(size(xi)), -1/4*ones(size(xi))]; 
        N2etaxi = [1/4*ones(size(xi)), -1/4*ones(size(xi)), 1/4*ones(size(xi)), -1/4*ones(size(xi))]; 
        
    elseif p == 2
        N    = [xi.*(xi-1).*eta.*(eta-1)/4, xi.*(xi+1).*eta.*(eta-1)/4, ...
            xi.*(xi+1).*eta.*(eta+1)/4, xi.*(xi-1).*eta.*(eta+1)/4, ...
            (1-xi.^2).*eta.*(eta-1)/2,  xi.*(xi+1).*(1-eta.^2)/2,   ...
            (1-xi.^2).*eta.*(eta+1)/2,  xi.*(xi-1).*(1-eta.^2)/2,   ...
            (1-xi.^2).*(1-eta.^2)];
        
        Nxi  = [(xi-1/2).*eta.*(eta-1)/2,   (xi+1/2).*eta.*(eta-1)/2, ...
            (xi+1/2).*eta.*(eta+1)/2,   (xi-1/2).*eta.*(eta+1)/2, ...
            -xi.*eta.*(eta-1),          (xi+1/2).*(1-eta.^2),   ...
            -xi.*eta.*(eta+1),          (xi-1/2).*(1-eta.^2),   ...
            -2*xi.*(1-eta.^2)];
        
        Neta = [xi.*(xi-1).*(eta-1/2)/2,    xi.*(xi+1).*(eta-1/2)/2, ...
            xi.*(xi+1).*(eta+1/2)/2,    xi.*(xi-1).*(eta+1/2)/2, ...
            (1-xi.^2).*(eta-1/2),       xi.*(xi+1).*(-eta),   ...
            (1-xi.^2).*(eta+1/2),       xi.*(xi-1).*(-eta),   ...
            (1-xi.^2).*(-2*eta)]; 
        N2xi  = [-1/2*eta.*(eta-1)/2,   1/2*eta.*(eta-1)/2, ...
                   1/2*eta.*(eta+1)/2,   -1/2*eta.*(eta+1)/2, ...
                   -eta.*(eta-1),         1/2*(1-eta.^2),   ...
                   -eta.*(eta+1),        -1/2*(1-eta.^2),   ...
                   -2*(1-eta.^2)];
               
        N2eta = [-xi.*(xi-1).*1/4,    -xi.*(xi+1)*1/4, ...
                    xi.*(xi+1).*1/4,     xi.*(xi-1)*1/4, ...
                   -(1-xi.^2)*1/2,      -xi.*(xi+1),   ...
                    (1-xi.^2)*1/2,      -xi.*(xi-1),   ...
                   -(1-xi.^2)*2];
               
        N2xieta = [(xi-1/2).*(2*eta-1)/2,   (xi+1/2).*(2*eta-1)/2, ...
                   (xi+1/2).*(2*eta+1)/2,   (xi-1/2).*(2*eta+1)/2, ...
                   -xi.*(2*eta-1),          (xi+1/2).*(-2*eta),   ...
                   -xi.*(2*eta+1),          (xi-1/2).*(-2*eta),   ...
                    4*xi.*eta];
               
        N2etaxi = [(2*xi-1).*(eta-1/2)/2,   (2*xi+1).*(eta-1/2)/2, ...
                   (2*xi+1).*(eta+1/2)/2,   (2*xi-1).*(eta+1/2)/2, ...
                   -2*xi.*(eta-1/2),        (2*xi+1).*(-eta),   ...
                   -2*xi.*(eta+1/2),        (2*xi-1).*(-eta),   ...
                    4*xi.*eta];
    else
        error('not available interpolation degree')
    end
    
elseif elem == 1
    if p == 1
        N = [1-xi-eta, xi, eta]; 
        Nxi = [-ones(size(xi)), ones(size(xi)), zeros(size(xi))];
        Neta = [-ones(size(xi)), zeros(size(xi)), ones(size(xi))];
        N2xi  = [zeros(size(xi)), zeros(size(xi)), zeros(size(xi))];
        N2eta = [zeros(size(eta)), zeros(size(eta)), zeros(size(eta))];
        N2xieta = [zeros(size(xi)), zeros(size(xi)), zeros(size(xi))];
        N2etaxi = [zeros(size(eta)), zeros(size(eta)), zeros(size(eta))];
        
    elseif p == 2
        kpa = 1-xi-eta;
        kpx = -1;
        N = [xi.*(2*xi-1), eta.*(2*eta-1),kpa.*(2*kpa-1), 4*xi.*eta, 4*eta.*kpa, 4*kpa.*xi ];
        Nxi = [4*xi-1, zeros(size(xi)), 4*kpa.*kpx-kpx, 4*eta, 4*eta.*kpx, 4*kpx*xi+4*kpa];
        Neta = [zeros(size(eta)), 4*eta-1, 4*kpa.*kpx-kpx, 4*xi, 4*kpa+4*eta*kpx,  4*kpx.*xi];
        N2xi = [4*ones(size(xi)), zeros(size(xi)), 4*ones(size(xi)), zeros(size(xi)), zeros(size(xi)), -8*ones(size(xi))];
        N2eta = [zeros(size(eta)), 4*ones(size(eta)), 4*ones(size(eta)), zeros(size(eta)), -8*ones(size(eta)), zeros(size(eta))];
        N2xieta = [zeros(size(xi)), zeros(size(xi)), 4*ones(size(xi)), 4*ones(size(xi)), -4*ones(size(xi)), -4*ones(size(xi))];
        N2etaxi = [zeros(size(eta)), zeros(size(eta)), 4*ones(size(eta)), 4*ones(size(eta)), -4*ones(size(eta)), -4*ones(size(eta))];
    else
        error('not available interpolation degree')
    end        
else
    error('not available element')
end
end
